function hfile = srf_RecalcNormals(hfile)
% SRF::RecalcNormals  - recalculate normals via triangles
%
% FORMAT:       [srf] = srf.RecalcNormals;
%
% No input fields
%
% Output fields:
%
%       srf         altered object

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin ~= 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);

% get coordinates and neighbors
crd = bc.VertexCoordinate;
tri = bc.TriangleVertex;
ncd = size(crd);
ntr = size(tri, 1);

% build edge vectors
csa = crd(tri(:, 2), :);
csb = crd(tri(:, 3), :) - csa;
csa = csa - crd(tri(:, 1), :);

% compute cross products (normals)
crs = cross(csa, csb, 2);

% compute length of normals
crl = sqrt(sum(crs .* crs, 2));

% scale to length 1 for averaging
crs = crs ./ crl(:, [1, 1, 1]);
crs(isnan(crs) | isinf(crs)) = 0;

% distribute over vectors
nrm = zeros(ncd);
for tc = 1:ntr
    for vc = 1:3
        nrm(tri(tc, vc), :) = nrm(tri(tc, vc), :) + crs(tc, :);
    end
end
crl = sqrt(sum(nrm .* nrm, 2));
nrm = nrm ./ crl(:, [1, 1, 1]);

% set in surface and store
bc.VertexNormal = nrm;
bvqxfile_setcont(hfile.L, bc);
